clear all;clc;close all;
load USPS.mat
%rename
training_digits = train_patterns;
training_labels = train_labels;
test_digits = test_patterns;
%%
%step 1-display the first 16 images in training digits
figure(1)
for i=1:16
 digit=reshape(training_digits(:,i),16,16); %reshape 
 subplot(4,4,i)
 imagesc(digit');
end
%%
%step 2-Print out the mean digits
training_averages=zeros(256,10);
figure (2)
for k=1:10
 onedigit=training_digits(:, training_labels(k,:)==1); %find all data of digit k
 training_averages(:,k)=mean(onedigit')'; %compute the mean of digit k
 reshapedigit=reshape(training_averages(:,k),16,16); %reshape 
 subplot(2,5,k)
 imagesc(reshapedigit');
end
%%
%step 3-the simplest classification computations
test_classification=zeros(10,4649);
%Calculate the distance between the test digit and the mean of each digit
for k=1:10
 test_classification(k,:)=sum((test_digits-repmat(training_averages(:,k),[1 4649])).^2);
end
test_classification_res=zeros(1,4649);
for j=1:4649
 [tmp, ind] = min(test_classification(:,j)); %find the index of mininum
 test_classification_res(j)=ind;
end
test_confusion=zeros(10,10);
%print out test_confusion
fprintf('test_confusion:\n')
for k=1:10
 tmp=test_classification_res(test_labels(k,:)==1);
 for j=1:10
 loc=find(tmp==j);
 test_confusion(k,j)=size(loc,2);
 fprintf('%5d',test_confusion(k,j));
 end
 fprintf('\n');
end

%%
%step-4 SVD-based classification computation.
left_singular_vectors=zeros(256,17,10);
for k=1:10
[left_singular_vectors(:,:,k),~,~] = svds(training_digits(:,training_labels(k,:)==1),17); %svd
end
test_svd17=zeros(17,4649,10);
for k=1:10
 test_svd17(:,:,k) = left_singular_vectors(:,:,k)' * test_digits; %compute coefficient
end
test_digits_rank_17_approximation=zeros(10,4649);
%compute approximation errors
for k = 1:10
 for j = 1:4649
 tmp = norm(test_digits(:,j) -left_singular_vectors(:,:,k)*test_svd17(:,j,k)); %calculate the approximation error
 test_digits_rank_17_approximation(k,j) = tmp;

 end
end
for j=1:4649
 [tmp,ind] = min(test_digits_rank_17_approximation(:,j)); %find the index of minimum
 svd_classification(1,j) = ind;
end
test_svd17_confusion=zeros(10,10);
%print test_svd17_confusion
fprintf('test_svd17_confusion:\n')
for k=1:10
 tmp=svd_classification(test_labels(k,:)==1);
 for j=1:10
 loc=find(tmp==j);
 test_svd17_confusion(k,j)=size(loc,2);
 fprintf('%5d',test_svd17_confusion(k,j));
 end
 fprintf('\n');
end